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ABSTRACT 

Eta Carinae is an ideal astrophysical laboratory for studying massive binary interactions and 
evolution, and stellar wind- wind collisions. Recent three-dimensional (3D) simulations set 
the stage for understanding the highly complex 3D flows in 77 Car. Observations of different 
broad high- and low-ionization forbidden emission lines provide an excellent tool to constrain 
the orientation of the system, the primary’s mass-loss rate, and the ionizing flux of the hot 
secondary. In this work we present the first steps towards generating synthetic observations 
to compare with available and future HST /STIS data. We present initial results from full 3D 
radiative transfer simulations of the interacting winds in rj Car. We use the Simplex algorithm 
to post-process the output from 3D SPH simulations and obtain the ionization fractions of hy- 
drogen and helium assuming three different mass-loss rates for the primary star. The resultant 
ionization maps of both species constrain the regions where the observed forbidden emission 
lines can form. Including collisional ionization is necessary to achieve a better description of 
the ionization states, especially in the areas shielded from the secondary’s radiation. We find 
that reducing the primary’s mass-loss rate increases the volume of ionized gas, creating larger 
areas where the forbidden emission lines can form. We conclude that post processing 3D SPH 
data with Simplex is a viable tool to create ionization maps for rj Car. 

Key words: radiative transfer - binaries: close - stars: individual: Eta Carinae - stars: mass- 
loss - stars: winds, outflows 


1 INTRODUCTION 

Eta Carinae (77 Car) is an extremely luminous (L Tota i > 5 X 10 6 L e ) 
colliding wind binary with a highly eccentric (e ~ 0.9), 5.54 year 
orbit (Davidson & Humphreys 1997; Damineli, Conti & Lopes 
1997; Hillier et al. 2001; Damineli et al. 2008a, b; Corcoran et al. 
2010). 7 / a , the primary of the system, is our closest (2.3 + 0.1 kpc, 
Smith 2006) example of a very massive star (~ 100 M 0 , Davidson 
& Humphreys 1997). A Luminous Blue Variable (LBV), p A has an 
extremely powerful stellar wind with M nA ^ 8.5 x 10 -4 M© yr _1 and 
Voo * 420 kms -1 (Hillier et al. 2001, 2006; Groh et al. 2012). Ob- 
servations over the last two decades indicate that t/a’s dense stellar 
wind interacts with the hotter, less luminous companion star t/ b and 
its much faster (Voo ~ 3000 kms -1 , Pittard & Corcoran 2002), but 
much lower density (M m ~ 10 -5 M 0 yr _1 ), wind (Damineli et al. 
2008a; Corcoran et al. 2010; Gull et al. 2009, 2011). These wind- 
wind interactions lead to various forms of time-variable emission 
and absorption seen across a wide range of wavelengths (Damineli 
et al. 2008a). 

Observational signatures that arise as a result of the wind- wind 
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interactions are important for studying 77 Car as they provide crucial 
information about the physical properties of the as-yet unseen p B 
and the system as a whole. Three-dimensional (3D) hydrodynami- 
cal simulations show that the fast wind of t/ b carves a low-density 
cavity out of the slower, denser inner wind of p A for most of the or- 
bit (Okazaki et al. 2008; Parkin et al. 2011; Madura & Groh 2012; 
Madura et al. 2012, 2013; Russell 2013). The same simulations in- 
dicate that the hot post- shock gas in the inner wind- wind interaction 
region (WWIR) gives rise to hard (up to 10 keV) X-ray emission 
that varies over the 5.54-year period. Together with the models, 
spatially unresolved X-ray (Hamaguchi et al. 2007; Henley et al. 
2008; Corcoran et al. 2010), optical (Damineli et al. 2008a, b), and 
near-infrared (Whitelock et al. 2004; Groh et al. 2010) observations 
have helped constrain the geometry and physical conditions within 
the inner WWIR. 

In addition to the ‘current’ interaction between the two winds 
that occurs in the inner regions (at spatial scales comparable to the 
semi-major axis length a ^ 15.4 AU « 0.0067 arc sec at 2.3 kpc), 
larger scale (« 3250 AU ^1.4 arcsec in diameter) 3D hydrody- 
namical simulations exhibit outer WWIRs that extend thousands 
of AU from the central stars (Madura et al. 2012, 2013, hereafter 
M12 and M13, respectively). Long-slit spectral observations of 
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77 Car with the Hubble Space Telescope/Space Telescope Imaging 
Spectrograph (HST/ STIS) reveal these spatially-extended WWIRs, 
seen via emission from multiple low- and high-ionization forbidden 
lines (Gull et al. 2009, 2011; Teodoro et al. 2013). 

Using a 3D dynamical model of the broad, extended [Fem] 
emission observed in 77 Car by the HST/ STIS, Ml 2 confirmed the 
orbital inclination and argument of periapsis that Okazaki et al. 
(2008) and Parkin et al. (2009) derived using X-ray data. More 
importantly, M12 broke the degeneracy inherent to models based 
solely on X-rays or other spatially-unresolved data and constrained, 
for the first time, the 3D orientation of 77 Car’s binary orbit. Ml 2 
find that the system has an argument of periapsis oj ^ 240° to 285°, 
with the orbital axis closely aligned with the Homunculus nebula’s 
polar axis at an inclination i « 130° to 145° and position angle on 
the sky PA « 302° to 327°, implying that apastron is on the ob- 
server’s side of the system and that tj b orbits clockwise on the sky. 
The dynamical model of M12 was based on 3D smoothed parti- 
cle hydrodynamics (SPH) simulations of 77 Car’s colliding winds. 
A simple radiative transfer (RT) code was used to integrate the op- 
tically thin [Fe 111 ] emission and generate synthetic slit-spectra for 
comparison to the available HST / STIS data. Although very suc- 
cessful, M 12 used a semi-analytic approach to compute the volume 
of wind material photoionized by r /^ . Furthermore, the fraction of 
Fe 2+ as a function of T was estimated assuming collisional ioniza- 
tion equilibrium and available ion fraction data. Due to the lack, at 
the time, of a suitable code, proper 3D RT simulations of 7 / b ’s ion- 
izing radiation were not performed. The location and strength of 
the [Fe in] emission was thus based on geometrical criteria, while 
in reality, the population of forbidden states depends on the local 
ionization state of the medium. 

The goal of this paper is to improve considerably the mod- 
eling approach of M12 by computing full 3D RT simulations 
of the effects of 77 B ’s ionizing radiation on 77 Car’s spatially- 
extended WWIRs. We accomplish this by applying the Simplex 
algorithm for 3D RT on an unstructured Delaunay grid (Ritzerveld 
& Icke 2006; Ritzerveld 2007; Paardekooper, Kruip & Icke 2010; 
Paardekooper et al. 2011; Kruip et al. 2010) to recent 3D SPH sim- 
ulations of 77 Car’s binary colliding winds (Ml 3). We use SimpleX 
to obtain detailed ionization fractions of hydrogen and helium at 
the resolution of the original SPH simulations. This should allow 
us to predict much more precisely where, and to what extent, var- 
ious observed forbidden emission lines form. This paper lays the 
foundation for future work aimed at generating synthetic spectral 
data cubes for comparison to data obtained with HST /STIS as part 
of a multi-cycle program to map changes in 77 Car’s extended wind 
structures across one binary cycle from 2009 through 2015 (Gull 
et al. 2011; Teodoro et al. 2013). Comparison of the observations 
to the models should ultimately lead to more accurate constraints 
on the orbital, stellar, and wind parameters of the 77 Car system, 
such as 77 a’s mass-loss rate and 77 B ’s temperature and luminosity 
(Mehner et al. 2010, 2012; M13). 

While we focus specifically on the case of 77 Car, the numer- 
ical methods in this paper can be applied to numerous other col- 
liding wind (e.g. WR 140, WR 137, WR 19, Fahed et al. 2011; 
Lefevre et al. 2005; Williams, Rauw & van der Hucht 2009) and 
dusty ‘pin wheel’ (WR 104, WR 98a, Tuthill, Monnier & Danchi 
1999; Monnier, Tuthill & Danchi 1999) binary systems. One of the 
biggest remaining mysteries is how dust can form and survive in 
such systems that contain a hot, luminous O star. Coupled with 3D 
hydrodynamical simulations, SimpleX simulations have the poten- 
tial to help determine the regions where dust can form and survive 
in these unique objects. 


In the following section we describe our numerical approach, 
including the SPH simulations, the SimpleX code, and the RT sim- 
ulations. Section 3 describes the results. A discussion of the results 
and their implications follows in Section 4. Section 5 summarizes 
our conclusions and outlines the direction for future work. 

2 CODES AND SIMULATIONS 
2.1 The 3D SPH Simulations 

The hydrodynamical simulations were performed with the same 
SPH code used in M13, to which we refer the reader for details. 
Optically thin radiative cooling is implemented using the Exact In- 
tegration scheme of Townsend (2009), with the radiative cooling 
function A (T) calculated using Cloudy 90.01 (Ferland et al. 1998) 
for an optically thin plasma with solar abundances. The pre-shock 
stellar winds and rapidly-cooling dense gas in the WWIR are as- 
sumed to be maintained at a floor temperature = 10 4 K due to pho- 
toionization heating by the stars (Parkin et al. 201 1). The same ini- 
tial wind temperature (T wind ) is assumed for both stars. The effect 
of T wind on the flow dynamics is negligible (Okazaki et al. 2008). 

Radiative forces are incorporated in the SPH code via an ‘anti- 
gravity’ formalism, the details of which can be found in M13 and 
Russell (2013). The individual stellar winds are parametrized us- 
ing the standard ‘beta- velocity law’ v(r) = Voo(l - R+/rf, where 
Voo is the wind terminal velocity, R+ the stellar radius, and (set 
= 1 ) a free parameter describing the steepness of the velocity 
law. Effects due to ‘radiative braking’ (Gayley, Owocki & Cran- 
mer 1997; Parkin et al. 2011), photospheric reflection (Owocki 
2007), and self-regulated shocks (in which ionizing X-rays from the 
WWIR inhibit the wind acceleration of one or both stars, leading 
to lower pre- shock velocities and lower shocked plasma tempera- 
tures, Parkin & Sim 2013), are not included. These effects are not 
expected to play a prominent role in 77 Car at the orbital phases near 
apastron considered in this work (Parkin et al. 2009, 2011; Russell 
2013; Ml 3). We include the more important velocity- altering ef- 
fects of ‘radiative inhibition’, in which one star’s radiation field 
reduces the net rate of acceleration of the opposing star’s wind 
(Stevens & Pollock 1994; Parkin et al. 2009, 2011). However, be- 
cause we fix the mass-loss rates in our anti-gravity approach, pos- 
sible changes to the mass-loss due to radiative inhibition are not 
included. These changes are not expected to be significant in 77 Car 
and should not greatly affect our results or conclusions (Ml 3). 

Using an xyz Cartesian coordinate system, the binary orbit is 
set in the xy plane, with the origin at the system centre-of-mass 
and the major axis along the x-axis. The two stars orbit counter- 
clockwise when viewed from the +z-axis. By convention, t = 0 
(0 = £/2024 = 0) is defined as periastron. Simulations are started 
at apastron and run for multiple consecutive orbits. Orbits are num- 
bered such that 0 = 1.5, 2.5 and 3.5 correspond to apastron at the 
end of the second, third, and fourth full orbits, respectively. 

The outer spherical simulation boundary is set at r - 100 a 
from the system centre-of-mass, where a = 15.45 AU is the 
length of the orbital semimajor axis. Particles crossing this bound- 
ary are removed from the simulations. The computational domain 
is comparable in size to past and planned HST /S TIS mapping ob- 
servations of the interacting stellar winds in 77 Car’s central core 
(~ +0.67" ^ +1540 AU, Gull et al. 2011; M12; Teodoro et al. 
2013). As demonstrated by Gull et al. (2011) and M12, 3D simu- 
lations at this scale are necessary for understanding and modelling 
the extended, time-variable forbidden line emission structures that 
are spatially and spectrally resolved by HST/ STIS. 
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Table 1 . Stellar, wind, and orbital parameters of the 3D SPH simulations 


Parameter 

7A 

7B 

Reference 

M* (Mo) 

90 

30 

H01; 008 

R* (Rq) 

60 

30 

H01; H06 

Twind (10 4 K) 

3.5 

3.5 

008; M13 

M (10 -4 M 0 yr" 1 ) 

8.5, 4.8, 2.4 

0.14 

G12; P09 

Voo (kms -1 ) 

420 

3000 

G12; P02 

P orb (days) 

2024 


D08a 

e 

0.9 


C01; P09 

fl(AU) 

15.45 


C01; 008 


Notes: M* and R+ are the stellar mass and radius. T w i n( j is the initial wind 
temperature. M and Voo are the stellar-wind mass-loss rate and terminal 
speed, respectively. P 0 rb is the period, e is the eccentricity, and a is the 
length of the orbital semimajor axis. 

References: C01 = Corcoran et al. (2001); HOI = Hillier et al. (2001); P02 
= Pittard & Corcoran (2002); H06 = Hillier et al. (2006); D08a = Damineli 
et al. (2008b); 008 = Okazaki et al. (2008); P09 = Parkin et al. (2009); G12 
= Groh et al. (2012). 

The total number of SPH particles used in the simulations is 
roughly between 5 x 10 5 and 9 x 10 5 , depending on the value of 
M nA . The adopted simulation parameters (Table 1) are consistent 
with those derived from the available observations, although there 
has been some debate on the exact present-day value of M r]A (see 
M13 for details). In an effort to better constrain t/ a ’s current M, 
M13 performed multiple 3D SPH simulations assuming different 
M rjA . We use the same naming convention as M13 when referring to 
the simulations in this paper for the different M r]A , namely, Case A 
{M T]a = 8.5 x 10 -4 M 0 yr -1 ), Case B (M^ = 4.8 x 10 -4 M 0 yr -1 ), 
and Case C (M T]A = 2.4 x 10 -4 M 0 yr -1 ). We discuss the effects of 
the three values of M T]A on the RT calculations in Section 3.2. 

2.2 The Simplex Algorithm for Radiative Transfer on an 
Unstructured Mesh 

The Simplex algorithm, conceived by Ritzerveld & Icke (2006), 
implemented by Ritzerveld (2007), and further improved by 
Paardekooper, Kruip & Icke (2010) and Kruip et al. (2010), is de- 
signed to solve the general equations of particle transport by ex- 
pressing them as a walk on a graph. At the basis of the method lies 
the unstructured grid on which the photons are transported. A given 
medium (e.g. a density or optical-depth field) is typically sampled 
with a Poisson point process and the resulting point distribution is 
used to tessellate space according to the Voronoi recipe: all points 
in a cell are closer to the nucleus of that cell than to any other 
nucleus. The Voronoi nuclei are then connected by a Delaunay tri- 
angulation. The grid is constructed to describe the properties of the 
underlying physical medium, through which the photons travel, in 
such a way that more grid points are placed in regions with a higher 
opacity. The result is a higher resolution in places where it is needed 
most, i.e. where the optical depth is highest. 

Photons are transported from node to node along the edges 
of the Delaunay triangulation, where each transition has a given 
probability. In one computational cycle, every nucleus in the grid 
transports its content to neighbouring nuclei, optionally absorbing 
or adding photons. Which neighbours are selected for transport 
depends on the specific process. Even though Simplex was orig- 
inally developed for application in cosmological radiative transfer, 
its properties are still well suited for our purpose. 

In Section 2.2.1 we present the construction procedure for the 


SimpleX RT mesh starting from the SPH particle distribution. Sec- 
tion 2 . 2.2 describes the processes that determine the ionization state 
of the gas, such as collisional- and photo-ionization and recombi- 
nation, plus the specifics of their implementation in SimpleX (for 
further details see Chapters 4 and 5 of Kruip 2011). 


2.2.1 Grid Construction 

The 3D SPH simulations provide the time-dependent 3D density, 
temperature, and velocity structure of 77 Car’s interacting winds on 
the spatial scales of interest, thus forming the basis of our model. 
As in Ml 2, the radiative transfer calculations are performed as post- 
processing of the 3D SPH output. The hydrodynamic influence of 
the radiation on the gas is thereby neglected. We do not expect this 
to have a very large influence on the results since the material pho- 
toionized by responds nearly instantaneously to its UV flux, i.e. 
the recombination time- scale is very small relative to the orbital 
time- scale, especially around apastron (Ml 2 ). 

The first step is to convert the SPH particle distribution to a 
SimpleX mesh. Since the density field is given by discrete particles, 
we might obtain an estimate of the density at any position in the 
domain using a typical SPH kernel function W(r,h) with 

P(r) = X m i w (\ r ~ rj\, h) (1) 

j 

where h is the smoothing length and is the mass of particle j. 
Using this kernel function one can then sample the data using the 
sampling functions described in e.g. Paardekooper, Kruip & Icke 
(2010) and Kruip et al. (2010). However, given the fact that the 
original data is already particle-based, it is more natural for us to 
use the SPH particles themselves as the generating nuclei for the 
Voronoi-Delaunay mesh. This leads to a more direct estimate of the 
density, given by the division of the particle mass by the Voronoi 
volume of its corresponding cell. Another advantage is that, due 
to pressure forces, the particles in an SPH simulation are in gen- 
eral positioned more regularly than for a pure Poisson process. Fi- 
nally, we note that with future applications of 3D time-dependent 
radiation-hydrodynamics in mind, a coupling of SimpleX with an 
SPH method is most natural when the radiation transport is applied 
directly to the SPH particles so that no spurious interpolation is 
needed. For these reasons we use every SPH particle as the nucleus 
of a Voronoi cell. This procedure yields density estimates that are 
less smooth than those obtained with typical kernel functions of the 
type of Equation (1), but guarantees mass conservation and repre- 
sents small scale structures in the density field more accurately. 

Figure 1 presents an example of the resulting SimpleX mesh 
and number density at apastron for a typical 3D SPH simulation 
of tj Car. The first row shows the original number density from 
the SPH simulation for slices in the xy-, xz-, and yz-planes for the 
Case A simulation. The SimpleX mesh (second row) reproduces 
everywhere the features present in the original SPH data. The re- 
sulting SimpleX number density (third row) follows extremely well 
the SPH one in shape, resolution, and value. 

Paardekooper, Kruip & Icke (2010) and Kruip et al. (2010) 
discussed how undersampling can have a negative effect on the 
outcome of an RT simulation that uses grid-based data. In princi- 
ple, the same problem holds for particle-based data in areas where 
sharp gradients are present in the number density of SPH particles. 
To ensure our results are not prone to such issues, we developed a 
method to increase the resolution of any sparsely sampled regions. 
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xy plane 


xz plane 


yz plane 



Log n [cm -3 ] 


1.8 5.8 9.8 13.7 


Figure 1 . Slices in the xy- (left column), xz- (middle column) and yz- (right column) planes through the 3D simulation volume for the Case A simulation at 
apastron. Rows show, from top to bottom, the original SPH number density distribution (log scale, cgs units), the Simplex mesh, and the resulting Simplex 
number density (same log scale, cgs units). The resolution of the Simplex mesh, as well the number density, follow well the resolution of the original SPH 
data. In the first column (i.e. the orbital plane) t/a is to the left and t/b is to the right. The length scales are shown under the top and bottom left panels. Note 
that the domain size in the Simplex snapshots is slightly smaller than that of the SPH simulations only because we have, for visualization purposes, removed 
the border points used to generate the Simplex mesh. 
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We begin with the triangulation of the SPH particles. For every 
tetrahedron in this triangulation that is larger than a tolerance vol- 
ume, an additional vertex is placed in the centre of the tetrahedron, 
and 1/5 of the mass of the four vertices that constitute the tetrahe- 
dron is given to the new vertex. This procedure is manifestly mass 
conserving and regularizes the mesh in low resolution regions. 

To study the influence this procedure has on the RT results, 
we performed test RT simulations with and without resolution en- 
hancements in normally sparsely- sampled regions. For these tests 
we used a characteristic 77 Car SPH simulation snapshot in which 
the densities span roughly ten orders of magnitude (i.e. the top row 
of Figure 1). If the RT results are sensitive to sharp gradients, we 
expect significant differences between the RT simulation with in- 
creased resolution and the original non-enhanced simulation. 

Although the difference is not large, the overall density in 
the affected regions is slightly enhanced by the new interpolation. 
However, we find that the overall shape and spatial extent of the 
ionized regions does not change. The ionization fraction in the lo- 
cal region most affected by the increased resolution is somewhat 
lower after the procedure though. This is consistent with the notion 
that the interpolation increases slightly the density in the higher- 
resolution region, resulting in a slightly higher recombination rate. 
The difference, however, is negligible. These results suggest that 
locally increasing the resolution does not change the overall shape 
and size of the ionized regions. We are therefore confident that our 
RT results are not susceptible to systematic effects related to strong 
gradients in the SPH particle number density (for further detail see 
Chapters 3 and 4 of Kruip 2011). 

2.2.2 Ionization State and Chemistry of the Gas 

In a gas with cross section for photoionization cr(v, v) at position x, 
the local photoionization rate, T P fx) (which gives the number of 
photoionizations per second per atom of species i in units [s -1 ]), is 
given by (Osterbrock & Ferland 2006): 

r° 47r/ v (f) 

r > i(x) = J (Ti(x , v) dv, (2) 

where J v (x) is the local mean intensity and the three species ca- 
pable of absorbing ionizing photons in the code are Hi, He 1 and 
He 11. For simplicity, we only consider the ionization of hydrogen 
and helium atoms. Implementing ionization processes in a numeri- 
cal code requires that the relevant equations be expressed in a dis- 
cretised form. In particular, we need to know the ionization rate per 
species in each cell of our computational grid. In Simplex, ioniz- 
ing radiation travels from cell to cell along the Delaunay edges. At 
the nucleus of each Voronoi cell, photons are taken away from the 
incoming radiation field and their energy is used to ionize the neu- 
tral atoms of that Voronoi cell. Given the number densities of these 
species (t7 H i, ^Hei and 7z He ii), and the path length through the cell /, 
the monochromatic optical depth of ionizing radiation r y is 

T v - (WhiO" HI + n HoI cr He I + ^Hell HellV (3) 

The total number of ionizations per unit time, Von for a cell with 
optical depth r y is then given by 

Von = f V r (v)[l - exp(-r y )] dv, (4) 

Jo 


where N y (v) is the number of ionizing photons per unit time stream- 
ing into the cell. To quantify how much of the resulting ionizations 
is due to a particular species, we use the contribution to the total 
optical depth of that species. The number of ionizations of species 
i per unit time is 

r°° 7" 

Mon, ;= Mon — dv (5) 

Jo T v 

Dividing by the number of neutral atoms of species i in the cell, V, 
then gives the spatially discretised equivalent of Equation (2) 


In numerical simulations involving radiation it is often nec- 
essary to approximate the continuous spectrum of radiation with a 
finite number of discrete frequency bins due to memory require- 
ments. The extreme (but often employed) limit of one single fre- 
quency bin is commonly referred to as the ‘grey approximation’. 
Although in the grey approximation all spectral information is lost, 
it is still possible to enforce the conservation of a quantity of impor- 
tance such as the number of ionizations per unit time or the energy 
deposition into the medium per unit time. For simplicity, we em- 
ploy the grey approximation in this work. 

The conservation of ionizations is accomplished by defining 
the effective cross section for species /, 07 as 

07, ; = 4x J dv /M ( 7 ) 

where N is the rate of ionizing photons per surface area defined by 

N = 4n f f dv. (8) 

Jo hv 

The photoionization rate is thus given by 

r P ,i = o- Li N. (9) 

In addition to photoionization, we include collisional ioniza- 
tion due to the interaction of free electrons and neutral atoms. As 
this is a kinetic process, the collisional ionization rate, T c , depends 
on the thermal state of the electrons and is given by 

r c = X r ‘ (r)n '> (i0) 

i 

where T^T) are the collisional ionization rates and n e is the electron 
number density. The total ionization rate is the sum of photo- and 
collisional ionization rates, T = T P + T c . 

The inverse process of ionization is recombination. This free- 
bound interaction of electrons and ions depends on temperature and 
number density of ions and electrons. The number of recombina- 
tions per unit time per hydrogen atom [s -1 ] is 

Ri - n^afT), ( 11 ) 

where afT) is the recombination coefficient of species i. 

We note that we use the ‘case B’ recombination coefficient 
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where the recombination transition to the ground- state is excluded 
under the assumption that the radiation associated with this transi- 
tion is absorbed nearby, resulting in a new ionization. This is re- 
ferred to as the ‘on-the-spot’ approximation. Details on the imple- 
mentation of these processes, as well as the various rates and cross 
sections used can be found in Chapter 5 of Kruip (201 1). 

Together, ionizations and recombinations determine the 
ionization- state of the gas, described by the following three cou- 
pled differential equations and three closure relations 


^hi = ^HII^HII - n m r m 

%el - ^Hen^Hell “ ^Hell^He! 

^Hell - ^Helll^HcIII ~ ^He IlI^He n 

rn = n m + nu II 

^He “ ^Hel + ^Hell + ^Helll 

n e = n m i + ft He ii + 2 n He iii- 


( 12 ) 


This set of equation does not have a general analytical solution 
and must be solved numerically. For this purpose we adopt the 
sub-cycling scheme described in Pawlik & Schaye (2008). In this 
scheme, ionizations and recombinations are evolved on a timescale 
that is smaller than the ionization or recombination timescales f ion 
and t rec . During a radiative transfer time-step, the ionizing flux is 
assumed to be constant, making the procedure manifestly photon- 
conserving. This allows for radiative time-steps A t n that are much 
larger than the dominant timescale governing the evolution of the 
ionization state. The sub-cycling time-step for both ionization and 
recombination is 


choice lies in the slow dynamical changes that the system under- 
goes around apastron. This ‘stable’ situation allows us to run RT 
simulations for a sufficiently long time without worrying about im- 
portant changes to the 3D structure of the system that occur around 
periastron (Okazaki et al. 2008; Parkin et al. 2011; M12; M13). 
Moreover, the HST/STIS mapping data currently in-hand to be 
modeled was taken at phases around apastron during 77 Car’s orbital 
cycle (Gull et al. 2011; Teodoro et al. 2013). Detailed modeling of 
future (late 2014 through early 2015) HST observations obtained 
across 77 Car’s periastron event is deferred to future work. 

We focus on the ionization of H and He due to r] B , assuming 
the same abundance by number of He relative to H as Hillier et al. 
(2001), 77 He / nu = 0.2. The reasons for this single-source approxi- 
mation are discussed in Section 2.3.1. We performed tests to deter- 
mine the correct time- step for accurate RT calculations of the ion- 
ization volumes and fractions, and find that a simulation time-step 
of ~ 3 s is required. The SPH output is post-processed with Sim- 
plex until the ionization state reaches an equilibrium value. This 
typically happens within ~ 3 months for the SPH snapshots inves- 
tigated. We thus set the total Simplex simulation time to 3 months. 
Because the gas is assumed to be initially almost fully neutral, this 
provides an upper limit on the time it takes before convergence is 
reached. Since this limit is well within the orbital time- scale around 
apastron, this is another indication that post-processing of the SPH 
simulations does not significantly alter our results. 


2.3.1 Influence of the Primary Star t/a 
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A?sub - ■ (13) 
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Because the procedure is analogous for each species, we give here 
only the explicit example for the integration step for hydrogen. At 
time t mb e (t n , t n + A t n ) the rate equation is given by 
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where the photoionization rate at / sub is 
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where T H and r are the photoionization rate and optical depth at the 
beginning of the sub-cycling and r (?sub) = rn^^/nm. By defining 
the photoionization rate in this way, the ionizing flux in the cell 
is constant during the radiative transfer time-step. This sub-cycling 
scheme becomes computationally expensive when Az sub «: A t n , but 
photoionization equilibrium is generally reached after a few sub- 
cycles. It is then no longer necessary to explicitly integrate the rate 
equation, but instead use the values of the preceding sub-cycle step. 
This way of sub-cycling ensures photon conservation even for large 
radiative transfer time-steps. 


2.3 Application of Simplex to 77 Car 

Since the Simplex calculations are performed as post-processing 
on the 3D SPH simulation output, we use snapshots correspond- 
ing to an orbital phase of apastron (Figure 1). The reason for this 


Detailed fitting of the optical and UV spectra of 77 Car by Hillier 
et al. (2001, 2006) and Groh et al. (2012) shows that for M r]A & 
8.5 x 10 -4 M© yr -1 , the region of fully ionized H around r/A extends 
radially ~ 120 AU, while the region of doubly-ionized He extends 
~ 0.7 AU, and that of singly-ionized He from ~ 0.7 to 3 AU. As- 
suming a constant spherical mass-loss rate, the density in the t] A 
wind is expected to fall off as r~ 2 . To explore the dependence of the 
position of the ionization fronts on the ionizing luminosity of t/a, 
we performed ID calculations using an equilibrium chemistry solu- 
tion where the ionization fractions of hydrogen and helium are set 
to their equilibrium values under the assumption that the incoming 
flux of ionizing photons is constant. As mentioned in Section 2.2.2, 
the ionization state of the gas is described by the first three equa- 
tions in set (12). We can derive the equilibrium equations by setting 
t/ri = ^Hei = %eii = 0- After some algebra this yields 


Xm 

= (1 +r H i/fl H n) 


•*HII 

- 1 - x m 


•*HeI 

= [1 + r He l/^HeII X (1 + r He II / ^He III)] 1 

(16) 

•*HeII 

- -^He ir He I / ^He II 


•XHelll 

- -^He II r He II / ^He III > 


where 

Xi is the fraction of species i and we have used n t 

= jq-TZj 


where j e (H, He). These equations are coupled by the free electron 
density given by the last equation in ( 12 ). 

Unfortunately, the set of equations (16) cannot be solved an- 
alytically due to the non-linear dependence on ionization fractions 
of the photoionization rate through the optical depth. More specif- 
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Figure 2. Neutral (x, solid lines) and ionized (1- x, dashed lines) fractions 
of H (log scale) as a function of distance for the ID RT simulations of p A 
described in Section 2.3.1 and luminosities of Top : 1.9 x 10 53 s -1 , Centre : 
1.91756xl0 53 s -1 , and Bottom: 1.91758xl0 53 s -1 . This behaviour signifies 
that the H ionization front is highly unstable and thus changes qualitatively 
with small perturbations in either the density or the luminosity. 


ically, the photoionization rate in a cell is given by Equation (6), 
where the monochromatic analog of Equation (4) is 

Mon = N y [l ~ exp(-r)], (17) 

Withr = (^HI^H^HI+^Hel^He^Hel + ^Hell^He^HellM Because of this 
non-linear dependence, the equilibrium fractions must be solved 
iteratively. If the iterative procedure converges, the neutral fractions 
are assigned to the cell under treatment and the flux is diminished 
by the number of absorptions during that time-step At (nmTmAt). 

The density profile used in the ID code is obtained from the 
expression p(r) = M^ A /4^r 2 v(r), where v(r) = Voo(l - R^/rf (see 
Ml 3, eqn. Al). The ID code simulates radiation traveling through 
spherically symmetric shells with a maximal radius of the simu- 
lation, « 1622 AU. The radiation is injected in the first shell and 
then travels outward until it is either absorbed or exits the last shell. 
For the results shown we used ~ 3 x 10 4 shells. We note that time- 
stepping is arbitrary because of the equilibrium chemistry. The only 
variable is therefore the luminosity of the source p A . 

Using M r]A ^ 8.5 x 10 -4 M©yr _1 , for luminosities below a 
critical value (1.9 x 10 53 s -1 ), the computational box is neutral and 
the Stromgren radii are confined to the central 0.7 AU (top panel, 


Figure 2). The H ionization front is located somewhere between 
the centre and the outside of the box for a very small range of lu- 
minosity values (centred around 1.91756 X 10 53 s -1 , centre panel, 
Figure 2). The slightest increase in luminosity results in a com- 
pletely ionized box, while further increases result only in a lower 
neutral fraction throughout the simulation volume (bottom panel, 
Figure 2). This behaviour is completely expected, however, for ion- 
ization fronts in power-law density profiles with powers less than 
-2/3 (Franco, Tenorio-Tagle & Bodenheimer 1990; Shapiro et al. 
2006). For such profiles, the circumstellar medium simply cannot 
support stable ionization fronts. 

For these reasons, constraining the ionization fronts in rj a’s 
wind to the values derived by Hillier et al. (2001, 2006) and Groh 
et al. (2012) using Simplex is practically impossible given the 
fronts’ intrinsically unstable nature. We realize this ID result is 
over-simplified since the instability is real in a pure H or H + He 
gas, but will disappear with the introduction of the myriad of spec- 
tral lines, mostly by Fe, that have a so-called ‘line blanketing’ ef- 
fect on the stellar spectra. However, the inclusion of such additional 
species and blanketing effects is beyond the scope of this paper. 

Given the above difficulties, the most sensible choice for an 
initial effort to model the ionized WWIRs is to omit t/ a ’s radia- 
tion altogether. This may seem an oversimplification at first, but 
there are several arguments for this approximation. First, for the 
high M T]a Case A and B simulations, ID CMFGEN models by Hillier 
et al. (2001, 2006) show that the primary source will sustain an ion- 
ized hydrogen region that spans roughly 240-260 AU in diameter, 
~ 0.04%-0.05% of the volume of the SPH and Simplex simula- 
tions. The same models show that the total sum volume of singly - 
plus doubly-ionized helium in the central rj A wind accounts for 
~ 10 _6 % and 10 _5 % of the Simplex simulation volume for Cases A 
and B, respectively. These volumes are too close to the primary 
and too small to directly affect the ionization fronts and fractions 
at the locations where the spatially-extended, high-velocity forbid- 
den line emission forms, especially at orbital phases near apastron 
(Gull et al. 2009, 2011; M12; M13). They may, however, influ- 
ence the ionization structure further away indirectly by reducing 
the opacity for photons from the secondary source. This may be 
especially true very close to periastron. We expect that this would 
primarily result in UV flux from t/ a penetrating the WWIR more 
easily, effectively increasing the ionized fraction on the far side of 
the primary source. For an observer on Earth though, this region is, 
at periastron, located behind t/ a and therefore likely obscured by 
the dense primary wind. 

Second, although extremely luminous, because it is en- 
shrouded by a dense, optically-thick wind, t/ a has a spectrum rep- 
resentative of a much cooler star than (Hillier et al. 2001, 2006). 
The effective temperature of t] A at optical depth r = 2/3 (r ~ 4 AU) 
is predicted to be ~ 9200 K for Cases A and B (Hillier et al. 2001, 
2006). The ionizing flux from r/ A is thus substantially diminished 
before reaching the WWIR, located ~ 20-22 AU from r/ A when the 
system is near apastron. Since essentially zero photons with ener- 
gies above 13.6 eY from p A reach the WWIRs on the apastron side 
of the system at times near apastron, omission of the p A source is 
a justifiable simplification when the focus is on forbidden emission 
lines with ionization potentials above 13.6 eV. 

One might try to argue that because the ionized hydrogen and 
helium volumes in the inner primary wind extend far enough to en- 
compass both stars and the WWIR at phases close to periastron, 
photons from t] A will also reach the apastron side of the simula- 
tion volume. This argument relies, however, on the assumption that 
the ionized regions are indeed spherical and therefore penetrate the 
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WWIR toward the secondary star. This assumption is likely incor- 
rect given the high density of the WWIR. In other words, we would 
be applying a model based on spherical symmetry to a region that 
clearly has a very asymmetrical geometry. 

The exception to the above arguments is the Case C simula- 
tion. In this instance, according to ID CMFGEN models (Hillier et al. 
2006; Ml 3), H is fully ionized in the pre-shock t/ a wind through- 
out the entire simulation domain. Moreover, the He n region in the 
inner 77 A wind extends radially ~ 120 AU. While neglecting the rf\ 
ionizing source in this case likely produces incorrect RT results in 
the regions of 77 A wind on the periastron side of the system, for our 
purposes, the situation is actually not so bad. First, we note that, 
based on previous works, the mass-loss rate of t] A is very likely not 
as low as the value assumed in the Case C situation (see arguments 
in e.g. Hillier et al. 2006; Parkin et al. 2009; Teodoro et al. 2012, 
2013; Russell 2013; M12; M13), and so we will not be using the 
Simplex results obtained here for Case C to model the observed 
broad, high-ionization forbidden line emission. Rather, in addition 
to investigating how a reduced M r]A affects the ionization structure 
on the apastron side of the system, we use Case C as an illustrative 
example to determine whether 77 B ’s ionizing radiation can penetrate 
the dense WWIRs and further affect the ionization state of 77^ s 
wind. Having H initially ionized in the pre- shock r/ A wind would 
primarily influence the ionization structure indirectly by reducing 
the opacity for photons from t] B (assuming they can penetrate the 
WWIR), increasing the ionized fraction of H in the pre- shock r/ A 
wind, but having little effect on the overall ionization volume. 

Regarding the He ionization structure in Case C, because the 
inner Hen region extends ~ 120 AU, the innermost WWIR pene- 
trates rj A s Hen zone, even at apastron. However, the total volume 
of this inner He 11 region is still only ~ 0.04% of the total Simplex 
simulation volume, again too small to directly affect the ionization 
fronts and fractions at the locations where the high-ionization for- 
bidden lines of interest form. Assuming He is neutral in the rj A wind 
at the start of the Simplex simulations also allows us to more easily 
determine whether He-ionizing photons from t] b can penetrate the 
WWIRs and affect the pre- shock rj A wind. If so, the primary effect 
will be an increased fraction of He 11 in the innermost r/ A wind, with 
little to no effect on the shape or extent of the He 11 ionization vol- 
ume. Moreover, since is thought to be an O- or WR-type star 
with r e ff - 36,000-41,000 K (Verner, Bruhweiler & Gull 2005; 
Hillier et al. 2006; Teodoro et al. 2008; Mehner et al. 2010), the 
number of photons it produces capable of ionizing He 11 to He 111 
is effectively zero. Thus, there will be no He 111 region created by 
77 B beyond the WWIR zone, even if 77 B ’s radiation can penetrate 
the dense post-shock gas. Therefore, even for simulation Case C, 
the neglecting of t] A s radiation is a justifiable simplification for the 
purposes of our work. The only caveat is that we do not account for 
any possible ionization of rj A s pre-shock wind to Hem by soft X- 
rays produced in the 420 km s -1 r} A shock. However, any such He 111 
region near the WWIR zone at times around apastron is very likely 
to be negligible in extent, if it exists at all, as evidenced by the ab- 
sence of any significant detectable He 11 T4686 emission in 77 Car 
during its spectroscopic high state (Mehner et al. 2011; Teodoro 
et al. 2012). 

Based on the above considerations, we neglect the t] A ionizing 
source in this work and focus on the influence of rj B . 



Figure 3. Slices in the xy orbital plane through the 3D simulation volume 
for the Case A simulation at apastron. The top plot shows the number den- 
sity (log scale, cgs units), while the bottom shows log temperature (K). The 
lower-density r/ B wind is shock-heated by the wind- wind interaction to tem- 
peratures up to ~ 10 8 - 5 K. On the periastron (left) side the denser t] A wind 
and compressed shells formed at periastron radiatively cool to T ~ 10 4 K. 


individual points randomly distributed about the sphere that defines 
the injection radius used in the SPH simulations for the wind of 77 B 
(30 R 0 ). We use a total of 50 source points, which is large enough 
to result in a nearly isotropic photoionizing source. We find that 
using more points has little effect on the RT results. The total lu- 
minosity is divided among all 50 points, forming the nodes of the 
grid that emit radiation. These nodes are also capable of absorbing 
any radiation emitted by neighbouring points in the SimpleX grid. 
Based on the work of Mehner et al. (2010); Verner, Bruhweiler & 
Gull (2005) and Ml 2, we assume for 77 B a total ionizing flux for hy- 
drogen and helium of 3.58 x 10 49 photons s -1 (3.02 x 10 49 capable 
of ionizing H 1 and 5.62 X 10 48 for ionizing He 1), consistent with an 
05 giant with T e ff ~ 40, 000 K (Martins, Schaerer & Hillier 2005). 


2.3.2 The Ionizing Source r] B 

For the RT calculations, we place a spherical ionizing source cen- 
tred at the location of 77 B . This ‘source’ is composed of a series of 


3 RESULTS 

To provide context for interpreting the RT results below, we begin 
with a brief description of the density and temperature structure of 
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the system in the orbital plane for the Case A simulation (Figure 3). 
Example number density slices in the xz andyz planes, plus number 
density slices in each plane for the Case B and C simulations, can 
be found in Figures 5, 7, and 9. Example slices showing tempera- 
ture for each case can be found in Ml 3. We focus on the Case A 
simulation as the M r]A assumed in this case most likely represents 
t/a’s current observed mass-loss rate (Groh et al. 2012; Ml 3). 

Pittard et al. (1998); Pittard & Corcoran (2002); Okazaki et al. 
(2008); Parkin et al. (2009, 2011); Madura et al. (2012); and Rus- 
sell (2013) showed that t/ b is between the observer and t/ a at apas- 
tron. Due to the highly eccentric binary orbit, t/ b spends most of 
its time near apastron (right side of panels in Figure 3), so that the 
relatively undisturbed wind of rj A is located on the far (periastron) 
side of the system. Across every periastron passage, the hot and 
low density wind of t/ b pushes outward into the slow, high-density 
t/a wind, leading to the formation of a thin, high-density wall sur- 
rounding the lower-density, trapped wind of t/ b (Parkin et al. 2011; 
Ml 3). This dense wall is accelerated to a velocity higher than the 
normal terminal velocity of t/ a ’s wind and expands creating a thin, 
high-density sheet of trapped primary wind material. During peri- 
astron passage the arms of the WWIR become extremely distorted 
by orbital motion as the binary stars move toward their apastron po- 
sitions. Moving back toward apastron, orbital speeds decrease and 
the t/ b wind cavity regains its axisymmetric conical shape (Okazaki 
et al. 2008; Parkin et al. 2011; Ml 3). 

Dense arcs and shells of tj a wind visible in the outer regions 
on the apastron side of the system in the top panel of Figure 3 high- 
light the fact that the binary has already undergone multiple orbits. 
Narrow cavities carved by t/ b in t/ a ’s dense wind during each peri- 
astron passage also exist on the periastron side of the system. Bor- 
dering these narrow cavities are the compressed, density-enhanced 
shells of primary wind formed as a result of the rapid wind- wind 
collision during each periastron. 

While the periastron side of the system is dominated by the 
dense wind of t/ a , the apastron side is dominated by the much 
lower-density, faster wind of t/ b , although arcs of compressed rj A 
wind also extend to the apastron side. These arcs are the remnants 
of the shells of t/a wind that flow in the apastron direction when 
t/ b is at periastron (Ml 3). The partially intact, most recent shell is 
visible just to the right of the centre of the image in the top panel. 

There is also a clear temperature asymmetry between the apas- 
tron and periastron sides of the system (bottom panel of Figure 3). 
The gas on the periastron side is much colder at T « 10 4 K. The var- 
ious wind-wind collisions on the apastron side produce large vol- 
umes of gas shock-heated to temperatures between 10 6 and 10 8 5 K. 
Because the gas on the apastron side is composed mostly of t/ b wind 
material of low density, it cools slowly and adiabatically, allowing 
it to remain hot throughout the 5.54-year orbital cycle. In contrast, 
the dense shells of post- shock primary wind cool radiatively very 
quickly down to T ~ 10 4 K (Parkin et al. 2011; Ml 3). The inner- 
most region of the system where the current WWIR is located, and 
the region where t/ b ’s wind collides with the latest ejected shell of 
primary wind, have the highest temperatures and are responsible for 
the observed time-variable X-ray emission (Hamaguchi et al. 2007; 
Okazaki et al. 2008; Corcoran et al. 2010; Parkin et al. 201 1). 

3.1 The Importance of Collisional Ionization 

The bottom panel of Figure 3 shows that the shocks induced by 
the violent wind- wind collisions heat the gas in the system to very 
high temperatures. Since the lower-density material from t/ b cools 
adiabatically, this gas remains extremely hot for most of the or- 


bit, at temperatures well above those where collisional ionizations 
become important (> 10 6 K). The collisional ionization fraction 
depends strongly on the temperature of the gas, which is in prin- 
ciple a function of the hydrodynamical motion, photo-heating, and 
multiple cooling terms. In this initial study, as a first approximation 
we use the temperature calculated by the SPH code to estimate the 
importance of collisional ionizations. In order to assess which pro- 
cess dominates, we performed a series of simulations with/without 
collisional-/photo-ionization. For brevity, we discuss here only the 
results for hydrogen for simulation Case A. Results for helium and 
Cases B and C are qualitatively similar. 

Figure 4 summarizes the results. The three panels represent 
the SimpleX output if we include, respectively, only collisional ion- 
ization, only photoionization, or both. For the only collisional ion- 
ization case, we assume collisional ionization equilibrium as an ini- 
tial condition to the RT, as described in Section 2.2.2. In this case, 
the overall ionization structure unsurprisingly follows the plot of 
the temperature in Figure 3. The cold, dense primary wind on the 
periastron side of the system, and the dense WWIRs of compressed 
primary wind that extend to the apastron side of the system (both 
in black in the first panel of Figure 4), remain mostly neutral. How- 
ever, hydrogen in the hot, lower-density regions of shocked sec- 
ondary wind are collisionally ionized (in blue and purple in the 
first panel of Figure 4). The t/ b wind in the outermost parts of the 
system is the most highly ionized due to the much lower density of 
the gas there, which results in less recombination. We also note in 
particular the two ‘fingers’ of highly ionized t/ b gas that extend into 
the primary wind, located at the bottom of the panel. More impor- 
tantly, we see that when only collisional ionizations are used, the 
dense WWIRs remain almost entirely neutral. 

In the case of only photoionizations from t/ b (middle panel of 
Figure 4), mainly the lower-density t/ b wind on the apastron side 
of the system is highly ionized (in yellow and orange). The lower- 
density hot fingers of t/ b wind trapped between the high-density 
walls of t/a wind show no ionization. These regions are effectively 
shielded from the ionizing flux of t/ b . Another important difference 
is the level of ionization in the t/ b wind. Photoionizations are ca- 
pable of reducing the fraction of neutral hydrogen by roughly four 
more orders of magnitude, compared to the case with only colli- 
sional ionization. Additionally, the t/ b wind closest to the centre of 
the simulation is the most highly ionized since, even though the 
density is higher there, the material is much closer to the luminous 
ionizing source. This is the exact opposite of what was observed 
in the case of only collisional ionizations. We also see that pho- 
tons from t/ b are capable of penetrating the innermost wall of t/ a 
wind material on the apastron side of the system, thus also highly 
ionizing it and the outer portions of t/ b ’s wind. Detailed examina- 
tion further shows that when photoionizations are used, the edges 
of the WWIRs facing t/ b can be significantly ionized (log f m < -3, 
Madura & Clementel 2014, in prep.). 

Using both collisional- and photo-ionizations results in a sit- 
uation that resembles a superposition of the first two panels (right 
panel of Figure 4). The t/ b wind on the apastron side remains highly 
ionized, but collisional ionization helps ionize the fingers of t/ b 
wind located at the bottom of the panel. Interestingly, the Hi in 
the fingers is slightly more ionized now compared to the case with 
only collisional ionizations. This is because, due to the now reduced 
opacity caused by including collisional ionization, photons from t/ b 
can more easily penetrate into the fingers and increase the over- 
all level of ionization. A similar effect is seen in/near the WWIRs, 
which are also slightly more ionized when both collisional- and 
photo-ionizations are included, compared to the case with only 
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Figure 4. Slices in the xy orbital plane through the 3D Simplex simulation volume for Case A at apastron showing the logarithm of the fraction of neutral 
hydrogen assuming Left : collisional ionization only, Centre : photoionization only, and Right : collisional- and photo-ionization. Including collisional ionization 
is necessary to ionize the small cavities in the primary wind and the ‘fingers’ of low-density 77B wind trapped between the higher-density walls of 77 A wind that 
form around periastron. 


photoionization. However, even when both collisional- and photo- 
ionizations are used, the dense t] A wind on the periastron side of the 
system remains neutral. 

Given all of these results, we consider both collisional- and 
photo-ionizations as necessary for any proper RT simulations of 
77 Car. The remainder of the results in this paper are based on sim- 
ulations that accordingly incorporate both. 

3.2 Overall Ionization Structure and Influence of M nA 

3.2.1 The orbital plane 

The top row of Figure 5 shows the Simplex number density in the 
orbital plane for simulation Cases A-C. As demonstrated by Ml 3, 
M T]a determines the overall shape of the WWIRs and the stability 
of the arcuate shells expanding on the apastron side of the system. 
Lowering M VA increases the opening angle of the shock cone cre- 
ated by t/b, increasing the volume of low-density t/ b wind. This 
is particularly noticeable in the size of the low-density fingers of 
77B wind that strongly reduce the volume of unperturbed primary 
wind. The lower the M^ A , the wider and more extended the fingers. 
In Cases B and C, the fingers extend to the back (periastron) side 
of t/a’s wind. The dense shells of r/ A wind on the apastron side are 
also more stable and remain intact longer for higher values of M VA 
(Ml 3). As a consequence, we expect that the 3D shape, position, 
intensity, and variability of the ionization depend strongly on M r]A . 

The middle and bottom rows of Figure 5 present, respec- 
tively, the computed fractions of H 1 and H 11 in the orbital plane. 
The WWIRs and high-density walls surrounding the lower-density 
trapped wind of tj b define the separation between the neutral and 
ionized-hydrogen regions. These high-density rj A wind structures 
are able to trap the hydrogen ionizing photons from t] B . We also see 
that as Mrj A decreases, the volume of ionized hydrogen increases 
greatly on both the apastron and periastron sides of the system. The 
larger fingers for Cases B and C allow the ionizing radiation from 


77B to penetrate into the low-density cavities that are carved within 
the back side of the primary wind every periastron passage. 

In the H 11 maps of Figure 5 it is possible to see a large fraction 
of neutral hydrogen at and to the periastron side of t] A . As described 
in Section 2.3.1, in reality, the hydrogen in this inner region should 
be ionized by t] A out to a radius of ~ 120 AU in Cases A and B, 
and everywhere in Case C. However, the absence of an t] A ionizing 
source in our simulations prevents this from occurring. Nonethe- 
less, the absence of an t] A source in our simulations reveals an im- 
portant result that may otherwise be missed, namely, that the high 
optical depth of the inner WWIR prevents any rj B ionizing photons 
from penetrating into the inner r/ A wind. The lack of any regions of 
ionized hydrogen in the unshocked primary wind on the periastron 
side of the system implies that regardless of the ionization structure 
of 77 a’s innermost wind, ionizing photons from 773 cannot penetrate 
the inner WWIR or significantly affect the dense r/ A wind on the 
periastron side of the system at times around apastron. 

Figure 6 illustrates the fractions of He 1, He 11 and He 111 in the 
orbital plane for the three M rjA simulations. Comparing to Figure 5, 
we see that the regions of He 111 correlate strongly with the regions 
of H 11. As expected, the regions of H 1 and He 1 are also correlated. 
The fully-ionized nature of helium in the lower-density 77B wind 
is due to the presence of large volumes of very-high-temperature 
shocked gas, plus the relatively close proximity of such gas to the 
hot, luminous 77B ionizing source. As with hydrogen, the helium 
in the denser primary wind is neutral. The trends as a function of 
M I]a seen in Figure 5 for the hydrogen ionization structure are also 
apparent in the plots of He 1 and He 111. This is a key result, as it 
implies that even with the lower M T]A of Case C, 77 B ’s He-ionizing 
radiation cannot penetrate significantly the dense WWIRs. 

The structure of He 11 (middle row of Figure 6) is more in- 
volved than that of He 1 and He 111. Interestingly, significant frac- 
tions of He 11 are principally located in the high-density walls of the 
WWIRs and outer edges of the dense fingers of r/ A wind that define 
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Figure 5. Slices in the xy orbital plane through the 3D Simplex simulation volume for the three different assumed M r/A . Columns illustrate, from left to right, 
Cases A to C. Rows show, from top to bottom, the Simplex number density (log scale, cgs units) and the computed fractions of Hi and Hu (log scale). 


the low-density fingers of rj B wind. Regions of lower-temperature 
unshocked t/ b wind also consist of mostly He n. The He n is seen 
primarily as a marker for the transition between the regions of He i 
and He hi. For this reason He n appears to be an excellent tracer for 
the high-density compressed post-shock r] A gas. 


3.2.2 The xz and yz planes 

To help the reader more fully appreciate the complex 3D structure 
of the simulation results, Figures 7-10 present slices showing the 
number density and H and He ionization structures in the xz and 
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Figure 6. Same as Figure 5, but with rows showing, from top to bottom, the computed fractions of Hei, Hen and Hem (log scale). In this and future plots of 
He m, the dashed circle marks the edge of the spherical computational domain. 


yz planes for each M r]A . The differences in ionization structure be- 
tween the three M TlA are even more apparent in these two planes. 
There is a clear left-right asymmetry in the density and ionization 
structure in each panel of the Figures. As in Figures 5 and 6, the 
regions of He hi correlate strongly with the regions of H n, while 
regions of He i correlate with those of H i. The overall volume of 
ionized material increases with decreasing M J]A . 

Figure 7 shows that the higher the value of M r]A , the smaller 
the wind cavities carved by t] b on the left (periastron) side of the 


system. They are practically invisible for Case A. As a result, hy- 
drogen and helium both appear neutral on the left in the Case A 
panels. Only the large t/ b wind cavity on the apastron side of the 
system is ionized in Case A. Figures 7 and 8 illustrate how the 
wind cavities on both the periastron and apastron sides of the sys- 
tem are much larger and remain hot for Cases B and C, resulting in 
well defined regions of ionized hydrogen and helium. 

The top row of Figure 7 also shows clear differences with M TlA 
in the density and fragmentation of the dense shell of r/ A wind ma- 
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Figure 7. Same as Figure 5, but for slices centered in the xz plane. 


terial on the right (apastron) side of the system. In Case A, the 
dense shell is more or less intact, while in Cases B and C it has 
fragmented considerably and started to mix with the lower-density 
77 B wind. This fragmenting shell produces an interesting hydrogen 
ionization structure on the right (apastron) side of the system that 
consists of an inner and outer region of low-density, highly-ionized 


77b wind (in yellow/orange, middle row of Figure 7) separated by 
a diffuse shell of denser, less-ionized r] A wind (in red). The middle 
row of Figure 8 again shows that the He n is located in the high- 
density walls of the WWIRs and outer edges of the dense fingers of 
77 a wind that define the low-density fingers of 77 B wind, thus tracing 
the compressed post-shock t] A gas. 
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Figure 8. Same as Figure 6, but for slices centered in the xz plane. 


Figure 9 shows that the cavities carved on the right (+y) side 
are always smaller than the ones carved on the left (-y), regardless 
of the value of M nA . However, the difference in cavity size between 
the +y and -y sides increases with decreasing M r]A . The larger cav- 
ities on the left (-y) side are also hotter, resulting in well-defined 
regions of ionized H and He concentrated on the left side. Finally, 
we see that the shells of dense, compressed rj A wind on the left are 
thicker and remain intact longer the higher the M^ A , reducing the 
overall volume of ionized material. We note that in reality, H should 
be ionized everywhere in Case C. 


4 DISCUSSION 

A major goal of this work was to improve upon the simple approach 
of M12 for computing the highly-ionized regions in the 77 Car sys- 
tem where various observed forbidden emission lines form. The 
ionization volumes in Ml 2 were based on geometrical criteria com- 
bined with a density threshold, and considered only photoioniza- 
tion of hydrogen due to t] B . Figure 1 1 shows an example of the 
photoionization region in the orbital plane for Case A at a phase 
near apastron, computed using the methods of Ml 2. The result is 
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Figure 9. Same as Figure 5, but for slices centered in the yz plane. 


a rather large Stromgren- sphere-like volume that predicts the dis- 
tance that H i ionizing photons from t/ b can travel. Comparing this 
to the Simplex results in Figures 5 and 6 we clearly see that the 
SimpleX method does a significantly improved job at computing 
the detailed structure of the various ionization volumes, including 
the penetration of t/ b ’s photons into the fingers of low-density wind 


carved within the optically-thick wind of rj A . The approach of Ml 2 
does not account for the extended WWIR arcs on the apastron side 
and thus overestimates the ionization extent in these regions. 

Effects due to collisional ionization and recombination were 
also not considered by M12. In addition to missing details in 
the ionization of the low-density fingers of tj b wind on the apas- 
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Figure 10. Same as Figure 6, but for slices centered in the yz plane. 


tron side, subtle variations in the ionization of t/aS wind and the 
WWIRs on the periastron side, due to recombination, are also ab- 
sent in the Ml 2 results. More importantly, the M 12 model does not 
compute any ionization fractions. The ion fraction is estimated us- 
ing tables and assuming collisional ionization equilibrium. In con- 
trast, the SimpleX method computes detailed ionization fractions 
for both hydrogen and helium. This provides estimates of the extent 
and magnitude of the ionization as a function of energy, previously 
unavailable information that is important for determining where the 


forbidden lines of different ionization potential form. Such infor- 
mation is also crucial for placing constraints on t/ b ’s ionizing flux. 

Ml 2 found that the observed broad forbidden line emission 
(Gull et al. 2009, 2011) depends strongly on M r]A and the ionizing 
flux from 7/ b . Ml 3 suggested that if the flux from t/ b remains con- 
stant, but M VA drops by a factor of 2 or more from an initial value 
of « 8.5 x 10 -4 M 0 yr -1 , then the photoionization region created by 
t/ b should increase considerably in size. The results of the SimpleX 
simulations in Figures 5-10 confirm this, implying that any recent 
decrease in M r]A (as speculated by e.g. Mehner et al. 2010, 2011, 
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Figure 11. Slice in the xy orbital plane showing the hydrogen photoioniza- 
tion region created by ?7 b (white = ionized) for the Case A simulation at 
apastron, computed using the approach in Ml 2. 


2012) should greatly change the spatial extent, location, and flux 
of the observed broad high-ionization forbidden emission lines. 
The results of this paper will be used in future work to compute 
synthetic slit-spectral observations of various forbidden lines (e.g. 
[Fen], [Fem]) for comparison to recent (Gull et al. 2011; Teodoro 
et al. 2013) and planned observations of 77 Car from HST/ STIS. 
Comparison of the synthetic and observational data can be used 
to place additional constraints on any recent changes in M r]A , im- 
portant for determining 7/ Car’s near- and long-term fates (Ml 3). 
The improved Simplex models will also be useful for refining the 
orbital orientation parameters obtained by Ml 2, and possibly also 
the stellar wind parameters and/or wind momentum ratio. 


5 SUMMARY AND CONCLUSION 

We showed that using Simplex for the post-processing of 3D SPH 
simulation output is a viable method to investigate the ionization 
state of the gas in a complicated colliding wind binary like 77 Car. 
Simplex provides detailed 3D results of the ionization volumes and 
fractions for various species of interest, in this case hydrogen and 
helium, and improves greatly upon simpler approaches such as that 
in Ml 2. Below we summarize our most important results. 

1 . The unstructured Simplex mesh reproduces everywhere the fea- 
tures present in the original 3D SPH simulation data, leading to 
a density map that is in excellent agreement with the original 
SPH one, even where sharp gradients are present. Simplex also 
preserves the high spatial resolution of the original SPH data. 

2. The inclusion of collisional ionization changes the ionization 


structure of hydrogen and helium most notably in the under- 
dense fingers of wind that form between the dense shells of 
77A wind created every periastron passage. Since these regions 
are typically shielded from 77 B ’s ionizing flux, including colli- 
sional ionization is important to achieve a more complete de- 
scription of the total ionized volume. 

3. Collisional ionization is important in reducing the total optical 
depth within regions composed of hot t] B wind that are heated 
to high temperatures by the various wind- wind collisions. This 
increases the efficiency of photoionization by 77B, allowing por- 
tions of the dense areas of post-shock t] A wind and WWIRs on 
the apastron side of the system to be ionized to varying degrees. 

4 . The SimpleX simulations show that the dense, innermost WWIR 
prevents the t] B ionizing radiation from penetrating far into the 
inner wind of t] A . At phases near apastron, hydrogen and helium 
ionization are concentrated on the apastron side of the system, 
with the periastron side consisting of mostly neutral t] A wind. 
However, as M riA is decreased, low-density fingers of ionized 77s 
wind penetrate the dense r] A wind on the periastron side. 

5. We find regions of He 111 correlate strongly with regions of Hu, 
while regions of H 1 strongly correlate with those of He 1. He 11 
is more complex and primarily marks the transition between the 
regions of He 1 and He 111. He 11 appears to be an excellent tracer 
for the dense, compressed post-shock rj A gas and WWIRs. 

6. Changing M nA results in quite different ionization volumes, with 
much more ionized gas present for lower M r]A . Significant vari- 
ations in ionization structure due to changes in M TlA are clearly 
apparent in the xz and yz planes as well as the orbital plane. 

7. The large apparent changes in ionization volume with decreas- 
ing M nA imply that any major decrease in M r]A should lead to 
significant observable changes in the spatial extent, location, 
and flux of the broad high-ionization forbidden emission lines. 
Future models based on the SimpleX results may be used to con- 
strain any such potential changes. 

In addition to helping us understand 77 Car’s recent mass-loss 
history, the past (Gull et al. 2011; Teodoro et al. 2013) and future 
ffiT/STIS spatial maps of 77 Car’s high-ionization forbidden emis- 
sion lines are a powerful tool that can potentially be used to better 
determine the nature of the unseen companion star 77s. Specifically, 
detailed 3D models of the forbidden line emission based on Sim- 
plex results like those presented here may allow us to place tighter 
constraints on 77 B ’s ionizing flux. This could then be compared to 
stellar models for a range of O (Martins, Schaerer & Hillier 2005) 
and WR (Crowther 2007) stars, providing a more accurate estimate 
of 77B ’s luminosity and temperature. 

While applied here to the specific case of 77 Car, SimpleX can 
be used to study numerous other colliding wind binaries or similar 
systems of astrophysical interest. Application of the SimpleX algo- 
rithm is also not limited to the post-processing of SPH simulation 
data, output from grid-based codes that use adaptive mesh refine- 
ment (AMR) may also be analyzed using SimpleX (Kruip 2011). 
And although SimpleX has been used in this paper to post-process 
hydrodynamical simulation data, this work helps set the stage for a 
future coupling of SimpleX with the SPH method in order to per- 
form 3D time-dependent radiation-hydrodynamics simulations of 
complex astrophysical phenomena (see e.g. Pelupessy et al. 2013). 
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